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We present an efficient quantum algorithm for simulating the evolution of a sparse Hamiltonian 
H for a given time t in terms of a procedure for computing the matrix entries of H. In particular, 
when H acts on n qubits, has at most a constant number of nonzero entries in each row/column, 
and \\H\\ is bounded by a constant, we may select any positive integer k such that the simulation 
\& , requires 0((log* n)t 1+1 ^ 2k ) accesses to matrix entries of H. We show that the temporal scaling 

cannot be significantly improved beyond this, because sublinear time scaling is not possible. 
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r\ 1 I. INTRODUCTION 

<d : 

There are three main applications of quantum computer algorithms: the hidden subgroup problem, with Shor's fac- 
torization algorithm one important example Q, search problems @, and simulation of quantum systems 0.L! Lloyd's 
method for simulating quantum systems 4] assumes a tensor product structure of smaller subsystems. Aharonov and 
Ta-Shma (ATS) p| consider the alternative case where there is no evident tensor product structure to the Hamiltonian, 
but it is sparse and there is an efficient method of calculating the nonzero entries in a given column of the Hamilto- 
, nian. Such respresentations of Hamiltonians can arise as encodings of computational problems, such as simulations 

■ of quantum walks 0, H, El EH or adiabatic quantum computations ^1 • 

Here we apply the higher-order integrators of Suzuki 0, ^| to reduce the temporal scaling from i 3 / 2 |j| or t 2 
to the slightly superlinear scaling ^ 1 + 1 / 2fc j where k is the order of the integrator and may be an arbitrarily large 
integer. We determine an upper bound on the number of exponentials required to approximate the evolution with 
a given accuracy. This enables us to estimate the optimal value of k, and therefore the fc-independent scaling in t. 
We than prove that, in the black-box setting, this scaling is close to optimal, because it is not possible to perform 
Oh, simulations sublinear in t. We also provide a superior method for decomposing the Hamiltonian into a sum for the 

■ problem considered by ATS, which dramatically reduces the scaling of n 2 |l4j or n 9 [j| to log* n for n qubits. This 
method is similar to "deterministic coin tossing" |Ie[. as well as Linial's graph coloring method |l6j| . 
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II. PROBLEMS AND RESULTS 



We commence with a statement of the problems that we consider in this paper and follow with the solutions that 
;_i " will be proven. 



Problem 1. The Hamiltonian is of the form H = Y^j=i^j- The problem is to simulate the evolution e %Ht by a 

sequence of exponentials e~ lHjt such that the maximum error in the final state, as quantified by the trace distance, 
does not exceed e. Specifically we wish to determine an upper bound on the number of exponentials, N cxp , required in 
this sequence. 

For this problem, the Hj should be of a form that permits e~ lHjt to be accurately and efficiently simulated for 
arbitrary evolution time t' . It is therefore reasonable to quantify the complexity of the calculation by the number 
of exponentials required. This problem includes the physically important case of simulating tensor product systems 
considered by Lloyd 4] , for which each Hj can be considered to be an interaction Hamiltonian. It also may be applied 
to the case where there is a procedure for calculating the nonzero elements in the columns 5]. In that case, each Hj 
is a 1-sparse Hamiltonian. The decomposition must be calculated, which requires additional steps in the algorithm. 

Our general result for Problem ^ is the following theorem. 

Theorem 1. When the permissible error is bounded by e, A cxp is bounded by 

N cxp < 2m5 2t (mr) 1+1 / 2t /e 1/2t , (1) 
for e < 1 < 2m5 k ~ 1 T, where r = and k is an arbitrary positive integer. 
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By taking k to be sufficiently large, it is possible to obtain scaling that is arbitrarily close to linear in r. However, 
for a given value of r, taking k to be too large will increase N cxp . To estimate the optimum value of k to take, we 
express Eq. as 



N cxp <2m 2 Te 2kln5+ln(mT/e)/2k 

The right-hand side has a minimum for 

k = round 



-ViogsCmr/e) + 1 

Here we have added 1 and rounded because k must take integer values. Adopting this value of k provides the upper 
bound 



N cxp < 4to 2 t e 2 V 1 " 51n (™^A) ; (2) 

for e < 1 < rnr/25. Eq. (J2J) is an expression for iVe X p that is independent of k. 

The scaling in Eq. J5J is close to linear for large tot. We show that this scaling is effectively optimal, because it is 
not possible to perform general simulations sublinear in r (see Sea llVfl . This result applies in the "black-box" setting, 
so it does not rule out the possibility that individual Hamiltonians have structure which allows them to be simulated 
more efficiently 

The second problem which we consider is that of sparse Hamiltonians. 

Problem 2. The Hamiltonian H has no more than d nonzero entries in each column, and there exists a black-box 
function f that gives these entries. The dimension of the space which H acts upon does not exceed 2™ . // the nonzero 
elements in column x are y\, . . . , yd' , where d! < d, then f{x, i) — (yi, H XtVi ) for i < d 1 , and f(x, i) — (x. 0) for i > d! . 
The problem is to simulate the evolution e~ tHt such that the maximum error in the final state, as quantified by the 
trace distance, does not exceed e. We wish to determine the scaling of the number of calls to f, Abb, required for this 
simulation. 

For each x, the order in which the corresponding yi are given can be arbitrary. The function / is an arbitrary 
black-box function, but we assume that there is a corresponding unitary Uf such that 

U f \x,i)\0) = \(f> x ,i)\yi,H Xyyi ), 

and we may perform calls to both Uf and lit. Here \4> x ,i) represents any additional states which are produced in the 
reversible calculation of /. 

ATS approached the problem by decomposing the Hamiltonian into a sum of Hj . We apply a similar approach to 
obtain the following theorem. 

Theorem 2. The number of black-box calls for given k is 

N hh e O ((log* n)d 2 b 2k {d 2 T) 1+1 ' 2k /e 1 ' 2k ) (3) 

with log* n = min{r| log^ n < 2} (the ^ indicating the iterated logarithm). 

The log* n scaling is a dramatic improvement over the n 9 scaling implicit in the method of ATS, as well as the n 2 
scaling of Childs [lj. 

III. HIGHER ORDER INTEGRATORS 

To prove Theorem ^ we apply the method of higher-order integrators. Following Suzuki [T^.[l3|. we define 

m 1 

5 2 (A)=n eHjA/2 n eHj ' A/2 > 

3=1 j'=m 

which is the basic Lie- Trotter product formula, and the recursion relation 

S 2k (X) = [S 2 k-2(PkX)} 2 S 2 k-2((l - 4pk)\)[S 2 k-2(Pk\)} 2 
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with p k = (4 - 4 1 /( 2 ' £ ~ 1) ) 1 for k > 1. Suzuki then proves that 12] 



exp [ J^HjX ] -& fc (A) 



eO(|A| 



(4) 



for |A| — ► 0. The parameter A corresponds to —it for Hamiltonian evolution. 

We first assess the higher-order integrator method in terms of all quantities t, m, k, and ||-ff||. Our result is 

Lemma 1. Using integrators of order k and dividing the time into r intervals, we have the bound 

< 5(2 x 5 k - 1 m,T) 2k+1 /r >k , 



exp ( -itJ2 H i ] - [S2k{-it/r)] 

3 = 1 



(5) 



for 



Amh^r/r < 1, 
(16/3)(2 x h k - x mrf k+1 /r 2k < 1. 



(6) 



Proof. Consider a Taylor expansion of both terms in the left-hand side (LHS) of Eq. Q. Those terms containing A 
to powers less than 2k + 1 must cancel because the correction term is 0(|A| 2fc+1 ), and terms with A 2fc+1 for k' > k 
must contain a product of 2k' + 1 of the Hj terms; thus 

(m \ oo L y 2fe' + l 

y: h,x = + £ A 2fc,+i y: c? n H 3l , 
j=l J k'=k 1=1 q=l 

The constants C\ and the number of terms depend on m and k. 

In order to bound and Ly , first consider the Taylor expansion of the exponential in the LHS of Eq. (0J. Because 
the operators Hj are in general noncommuting, expanding (Hi + • • • + H m ) 2k +1 yields m 2k +1 terms. Therefore the 
Taylor expansion contains m 2k +1 terms with A 2fc +1 . These terms have multiplying factors of l/(2fc' + 1)! because 
this is the multiplying factor given by the Taylor expansion of the exponential. 

To place a bound on the number of terms in the Taylor expansion of S2fe(A), note that S2k{X) consists of a product 

of 

2(m- l)5 fe_1 + 1 

exponentials. The expansion for S2k{t) may be obtained by expanding each of the exponentials individually. There 
will be no more than 

[2(m-l)5 fe - 1 + l] 2fe ' +1 

terms with A 2fe +1 . Because \pk\ < 1 and |1 — 4pk\ < 1, the multiplying factors for each of these terms must be less 
than 1. 

Each Hj satisfies \\Hj\\ < \\H\\ Defining A = \\H\\, and using standard inequalities we obtain 



oo L k , 2k'+l 

e^ ,+i e^' n 

k' = k 1=1 q=l 



< E |AA| 2fc ' +1 L fc < 

k'=k 

< E |AA| 2fe ' +1 {m 2fc ' +1 + [2(m - l)b k ~ l + l] 2fe ' +1 } 



k'=k 
oo 



< 2 E |AA| 2fc ' +1 [2m5 



fc _ ll2fc , +1 _ 2|2m5 fc - 1 AA| 2fc + 1 



1 - \2m5 k - 1 XA\ 2 



k'=k 



Therefore we obtain the inequality 



CXp [ XJ2 H 3 ] -S2k(\) 



< (8/3)|2m5 fc - i AA| 2,c+1 , 
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provided \2m5 k 1 AA| < 1/2. Substituting A = —it/r where r is an integer, and taking the power of r, gives the error 
bound 



exp 




[S 2k (-it/r)y 



< [l + (8/3)(2m5* _1 Ai/r) 



for Am&^kt/r < 1. This may alternatively be expressed as in Lemma ^ 



(7) 



□ 



By placing limits on the norm of the difference in the unitaries, we limit the trace distance of the output states. 
This is because 

\\Ux-UtW > \\Ui\il>)-U 2 \il>)\\ 



> -Tr 
~ 2 



uM){i>\ul-uM{i>\ul 

= d (wk^.wx^) , 

with D the trace distance. We now use this to prove Theorem ^ 
Proof. ( of Theorem QJ) Let us take 



= [4 



x o 



fc-l/2 



(mr) 



l + l/2fc / £ l/2fen 



(8) 



Given the restrictions e < 1 < 2m5 fc_1 r, it is easily seen that Eqs. © hold. In addition, the right-hand side of Eq. (J5J 
does not exceed e, so the error can not exceed e. 

Because the number of exponentials in S2&(A) does not exceed 2m5 
in Eq. ©, then we find that 



1 we have N exp < 2m5 k 1 r. If we take r as 



AW < 2m5 2fc (mr) 1+1 / 2fe /e 1/2fc . 



(9) 



Here the multiplying factor has been changed to take into account the ceiling function. Hence the order scaling is as 
in Eq. Q. □ 

This result may be used for any case where the Hamiltonian is a sum of terms that may be simulated efficiently. It 
may therefore be applied to the case of tensor product systems, where the individual Hj are interaction Hamiltonians. 
It can be also used for cases of the type of Problem 2, where the Hamiltonian is sparse. In this case we have the 
additional task of decomposing the Hamiltonian into a sum. 



IV. LINEAR LIMIT ON SIMULATION TIME 



We have shown that the simulation of any Hamiltonian may be performed arbitrarily close to linearly in the scaled 
time t. We now show that the scaling cannot be sublinear in r, provided the number of qubits can grow at least 
logarithmically with respect to r. The result is 

Theorem 3. For all positive integers N there exists a row- computable 2-sparse Hamiltonian H such that simulating 
the evolution of H for scaled time r = irN/2 within precision 1/4 requires at least t/2-k queries to H. 

Here a row-computable Hamiltonian means one where there is a method for efficiently calculating the nonzero 
elements in each row. 

Proof. The idea is to construct a 2-sparse Hamiltonian such that the simulation of this Hamiltonian determines the 
parity of A^ bits. It has been shown that the parity of A^ bits requires A^/2 queries to compute within error 1/4 [T^ |: 
therefore the Hamiltonian can not be simulated any more efficiently. 

First consider a Hamiltonian H acting on orthogonal basis states |0), . . . , \N), for which the nonzero matrix entries 
are 

(j + l\H\j) = (j\H\j + 1) = ^(N~j)(j + 1)/2. 

This Hamiltonian is equivalent to a J x operator for a spin A^/2 system, with the |j) being J z eigenstates. It is therefore 
clear that e~ mH \Q) = \N) and = N/2. 
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FIG. 1: Graph representing example the Hamiltonian in the proof of Theorem|3| States are represented by ellipses, and nonzero 
elements of the Hamiltonian are indicated by lines. The sequence of states \kj,j) with ko = is indicated by the solid line. 



Now we construct an augmented version of the above Hamiltonian, that corresponds to a graph with two disjoint 
lines with weights as above, where the lines "cross over" at the positions where bits X%, . . . ,Xn are 1. We add an 
ancilla qubit so the Hamiltonian H acts on basis states |0, 0), . . . , |0, N), \1, 0), . . . , |1, N). The nonzero matrix entries 
of H are 

{k',j + l\H\k,j) = (k, j\H\k',j + 1) = y/(N + l)/2 

for values of k and k' such that k © k! = Xj + i (where © is XOR). 

Thus, if Xj+i is zero, then there is a nonzero matrix element between \0,j) and |0, j + 1), as well as between \l,j) 
and |l,j + 1). If Xj+i is equal to 1, then the nonzero matrix elements are between |0,j) and + 1), as well as 
|l,i) and |0, j + 1). We may determine a sequence of bits ko, . . . , k^ such that kj © fcj+i = Xj + \. The Hamiltonian 
acting on the set of states \kj,j) will then be identical to that acting on the states \j) with the original Hamiltonian. 
It is therefore clear that e'^ H \k o ,0) = \k N ,N). 

The graph corresponding to a Hamiltonian of this type is shown in Fig. ^ The system separates into two distinct 
sets of states which are not connected. If the system starts in one of the states on the path indicated by the solid 
line, it can not evolve under the Hamiltonian to a state on the dotted line. From the definition of the kj, if fco = 0, 
then kj is the parity of bits X\ to Xj, and in particular fc^v gives the parity of all N bits. Thus if we start with the 
initial state |0, 0) and simulate the evolution e~ l7rH , we obtain the state | fcjv, N), where k^ is the parity of the N bits 
Xi, . . . , Xjy. Thus measuring the state of the ancilla qubit will give the parity. 

Let us denote the final state obtained by the simulation by \ip), and the reduced density operator for the ancilla by 
p anc . If the error probability is no less than 1/4, then D(p anc , |fc/v)(fc/v|) > 1/4, which implies that 

D(\i,){il;\,\k N ,N){k N ,N\) > 1/4. 

Hence, if there are no more than N/2 queries to the Xj, the error in the simulation as quantified by the trace distance 
must be at least 1/4. 

Each query to a column of H requires no more than two queries to the Xj (for column j we require a query to Xj 
and Xj + i). Thus, if there are no more than N/4 queries to H, then there are no more than N/2 queries to the Xj. 
In addition, the scaled time for the simulation is 

r = ||iT||t = txN/2. 

Thus the simulation of H requires at least N/4 — t/2tt queries to obtain trace distance error no more than 1/4. □ 

The form of this result differs slightly from that in the previous section, in that the cost is specified in terms of 
the number of queries to the Hamiltonian, rather than the number of exponentials. It is straightforward to show the 
following result for the number of exponentials. 

Corollary 1. There is no general integrator for Hamiltonians of the form H = Hi + Hi such that (trace distance) 
error < 1/4 may be achieved with the number of exponentials N cxp < t/2tt. 

By general integrator we mean an integrator that depends only on r, and not the Hamiltonian. 

Proof. We take H as in the preceding proof. This Hamiltonian may be expressed in the form H = Hi + H<i by 
taking H\ to be the Hamiltonian with (k',j + \\H\\k,j) nonzero only for even j, and Hi to be the Hamiltonian with 
(k',j + l\Hi\k,j) nonzero only for odd j. Each query to the H^ requires only one query to the Xj. For H\ {Hi), 
determining the nonzero clement in column j requires determining which of j and j + 1 is odd (even) , and performing 
a query to the corresponding Xj or Xj + \. 

Both Hi and Hi are 1-sparse, and therefore may be efficiently simulated with only two queries [ItI Hsj . If -/V cxp < 
t/2tt, the total number of queries to the Xj is no more than t/tt. Taking t — it and \\H\\ = N/2, the number of 
queries is no more than N/2. However, from Ref. |19j the error rate can be no less than 1/4. Hence error rate < 1/4 
can not be achieved with 7V oxp < t/2tt. □ 
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V. EFFICIENT DECOMPOSITION OF HAMILTONIAN 



Next we consider the problem of simulating general sparse Hamiltonians, as in Problem|2] Given that the dimension 
of the space does not exceed 2™, we may represent the state of the system on n qubits, and x and y may be n-bit 
integers. The real and imaginary parts of the matrix elements will be represented by n' bit integers (for a total of 2n' 
bits for each matrix element), where ml must be chosen large enough to achieve the desired accuracy. 

In order to simulate the Hamiltonian, we decompose it into the form H = Ylj=i-Hj, where each Hj is 1-sparse 
(i.e., has at most one nonzero entry in each row/column) If Hj is 1-sparse then it is possible to directly simulate 
exp(-iHjt) with just two black-box queries to Hj 0, 0. Since the value of m directly impacts the total cost of 
simulating H, it is desirable to make m as small as possible. The size of the sum may be limited as in the following 
lemma. 

Lemma 2. There exists a decomposition H = J^Li-Hj, where each Hj is 1-sparse, such that m = 6d 2 and each 
query to any Hj can be simulated by making 0(log* n) queries to H . 

Proof. From the black-box function for H , we wish to determine black-box functions for each Hj that give the nonzero 
row number, y, and matrix element corresponding to each column x. This black-box for Hj is represented by the 
function g(x,j), with output (y, (Hj) x<y ). If there is no nonzero element in column x, the output is (x,0). 

Intuitively, it is helpful to consider the graph Gh associated with H whose vertex set is {0, 1}™. Each vertex 
corresponds to a row or column number, and there is an edge between vertex x and y if the matrix element H XyV is 
nonzero. As H is Hermitian we take the graph to be undirected. We wish to determine an "edge-coloring" of G#, 
which is a labeling of the edges such that incident edges have different colors. Each edge color, j, then corresponds 
to a different Hamiltonian Hj in the decomposition of H. 

The basic idea is as in the following labeling scheme, where the labels are indexed from the set {1, ... , d} 2 . We take 
f y to be the y-component of /; then f y (x, i) gives the i th neighbor of vertex x in the graph. Let (x, y) be an edge of 
Gh such that y = f y (x,i) and x = f y (y,j)- Thus edge (x,y) is labeled with the ordered pair for x < y, or (J,i) 
for x > y. This labeling is not quite an edge-coloring; for w < x < y it is possible for edges (w, x) and (x, y) to both 
have the label That will be the case if y and w are the i th and j th neighbors of x, respectively, and x is the i th 

neighbor of w and the j th neighbor of y. To ensure that the labels are unique, we add the additional parameter v, so 
the label is (i,j, u). 

We assign v via a method similar to "deterministic coin tossing" [l5j . We set Xq = x, then determine a sequence 
of vertices 

4 0) <4 0) <4 0) <--. 

such that Xi?x = fy(x[°\i) and f y (xfvi,j) — x[°\ That is, the edges (x[°\ xj?^) are labeled (i,j,v), with the same 
values of i and j for each edge. We need to choose values of v for the edges such that the same value is never repeated 
in this chain. 

A typical chain may have only two elements; however, there exist Hamiltonians such that long chains may be 
formed. In the case that the chain is long, we do not determine it any further than x^ +1 . Here z n is the number of 
times we must iterate I i— > 2[log 2 T\ (starting at 2") to obtain 6 or less. This quantity is of order log* n, and for any 
realistic problem size z n itself will be no more than 6 [2l|. 

Now we determine a second sequence of values x\ j . This sequence is taken to have the same length as the first 
sequence. For each x[ ^ and ij? l5 we determine the first bit position where these two numbers differ, and record 
the value of this bit for x\ , followed by the binary representation of this position, as x^\ The bit positions are 
numbered from zero; that is, the first bit is numbered 00 ... 0. If xf^ is at the end of the sequence, we simply take 
x9^ to be the first bit of x, , followed by the binary representation of 0. There are 2™ different possible values for 
each of the x[°\ and 2n different possible values for each of the xf . 

From the definition, each xf 3 ' is unique. Also x^ must differ from This is because, even if the positions of 

the first bit where xf^ differs from x^ and x^ differs from x|9 2 are identical, the value of this bit for x^ will of 
course be different from the value for x£\. As the xf^ contain both the position and the value of the bit, x^ must 
differ from 5c|i\. 

There is a subtlety when xj^ is at the end of the sequence. Then x^ contains the first bit of x|?^, and the 
position of the first bit which differs is taken to be 1. In that case, if x^ differs from x^ at the first bit (so the bit 
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TABLE I: Example values of x\ p under our scheme for calculating v. The value of v obtained is in the upper right, and is 
shown in bold. For this example n — 18 and z n = 4. The values in italics are those that may differ from (there are no 

corresponding values for the bottom row). 



l\p 





1 


2 


3 


4 





001011100110011010 


000001 


0100 


000 


000 


1 


010110101010011011 


000010 


1100 


100 


100 


2 


011011101110101101 


000000 


0001 


000 


000 


3 


101011101011110100 


010001 


1001 


100 


100 


4 


101011101011110101 


000001 


0000 


000 


000 


5 


111000010110011010 


100000 


1000 


100 


100 



positions recorded in x^ and are identical), then the bit values which are recorded in x\ and x^_ 1 must be 
different. Thus it is still not possible for x^ to be equal to 

We repeat this process until we determine the sequence of values xf n \ We determine the from the Xj in 

exactly the same way as above. At each step, x[ p ' differs from xfi x for exactly the same reasons as for p = 1. As we 

go from p to p + 1, the number of possible values for the x| p ' is reduced via the mapping k > 2|~log 2 k~\ . Due to our 

(*») 



choice of z n , there are six possible values for x ' 



Now if w < x with x — f y (w, i) and w — f y (x,j), then we may set = w and perform the calculation in exactly 
the same way as for x in order to determine wi Zn \ If the chain of a^ ' ends before z n , then the x\ p ^ will be the same 
as the tw^Lj. In particular x Zn ^ will be equal to w[ Zn \ so it is clear that will differ from x Q Zn \ 

On the other hand, if there is a full chain of x^ up to x z + i, then the chain for w will end at w z +1 , which is 
equivalent to xf} ■ Then w z +1 will be calculated in a different way to x^J , and may differ. However, w^J will be 
equal to x^_ v At step p, w^j_ p+1 will be equal to x^_ p . In particular, at the last step, w[ z '^ will be equal to x Q Zn . 
Thus we find that w Zn •* again differs from Xg Z " . 

(z ) (z ) 

As x ™ has this useful property, we assign the edge (x,y) the color v), where v = x " . Due to the properties 
of the above scheme, if the edge (w,x) has the same values of i and j as (x,y), it must have a different value of v. 
Therefore, via this scheme, adjacent edges must have different colors. 

Now we describe how to calculate the black-box function g using this approach. We replace j with to 
reflect the labeling scheme, so the individual Hamiltonians are Huj v y The black-box function we wish to calculate 
is g(x,i,j, v). We also define the function T(x,i,j) to be equal to the index v as calculated in the above way. There 
are three main cases where we give a nontrivial output: 

1. f y (x,i) = x, i = j and v = 0, 

2. f y (x,i) > x, fy(f y (x,i),j) = x and T(x,i,j) = v, 

3- f y {x,j) < x, fy(f y (x,j),i) = x and T(f y (x, j), i, j) = v. 

In cases 1 and 2 we return g[x,i,j,v) = f(x,i), and for case 3 we return g(jx,i,j,v) = f(x,j); in all other cases we 
return g (x, v) = {x, 0). 

Case 1 corresponds to diagonal elements of the Hamiltonian. We only return a nonzero result for v = 0, in order to 
prevent this element being repeated in different Hamiltonians H^ j V y Case 2 corresponds to there existing a y > x 
such that y is the i th neighbor of x and x is the j th neighbor of y. Similarly Case 3 corresponds to there existing a 
w < x such that w is the j th neighbor of x and x is the i th neighbor of w. The uniqueness of the labeling ensures 
that cases 2 and 3 are mutually exclusive. 

As there are d possible values for i and j, and v may take six values, there are 6d 2 colors. Thus we may take 
?7i = 6d 2 . In determining v, we need a maximum of 2(z n + 2) queries to the black-box; this is of order log* n. □ 

To illustrate the method for determining v, an example is given in TableUfor {x[ p) } and Table|TT|for {w[ p} }. Fig.H 
shows a portion of the graph corresponding to the values in Tables III and ITU In the Tables n = 18, so there are 2 18 
possible values in the first column. Then there are 36 possible values in the second column, 12 in the third, 8 in the 
fourth and 6 in the fifth. Thus z n is equal to 4 in this case, and the sequence of x^ is determined up to x^ . 
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TABLE II: Example values of w\ under our scheme for calculating v. The value of v obtained is in the upper right, and is 
shown in bold. For this example n — 18 and z n = 4. The values in italics are those which may differ from x^} l . 



l\p 





1 


2 


3 


4 





000010010110111001 


000010 


1100 


100 


100 


1 


001011100110011010 


000001 


0100 


000 


000 


2 


010110101010011011 
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FIG. 2: A portion of the graph for the example given in Tables IJ and ITTI The vertices w, x, y, etc each have i = 1 and j = 3 for 
the edge labels, so it is necessary for the v to differ to ensure that adjoining edges have distinct labels. The xf and wf^ which 
the vertices correspond to are also given. The numbers in the first columns of Tables HI and ITTI are the binary representations of 
the vertex numbers given here. 



In Table [TJ the values of Xj are given, and the elements in the first column are values of xf^ . As an example of 
calculation of x\ , note that Xq "* differs from x^ in the second bit position. The second bit for Xq '' is 0, so this is 
the first bit for Xq . We subtract 1 from the bit position to obtain 1, and take the remaining bits of x^ to be the 
binary representation of 1. For the case of x!; '', this is the end of the chain, so we simply take to be the first bit 
of which is 1, and the binary representation of 0. 

In Table ITTI the values of are given, where these are calculated from a w < x such that x = f y (w,i) and 

w = f y (x,j). The example given illustrates the case where the sequence of wj® (with — x,_,) ends before the 

sequence of x, . In this case, we find that the differences propagate towards the left, but we still have Xq Z "' ) = wf . 

Thus different values of v are obtained, as expected. For x we obtain v = X(p = 000, and for w we obtain v — — 
100. 

We can use Lemma [21 to prove Theorem [2] 

Proof, (of Theorem^ Overall the number of Hamiltonians H^j^ in the decomposition is m — 6gP. To calculate 
g(x, i,j, v), it is necessary to call the black-box 2(z n + 2) times. 

To simulate evolution under the Hamiltonian Haj tV \, we require g to be implemented by a unitary operator U g 
satisfying 

U g \x,i,j,is)\0) = \x,i,j,v)\y,{H, hJtV ) x .y). 

As discussed above, the function / may be represented by a unitary Uf, using this unitary it is straightforward to 
obtain a unitary U g such that 

U g \x,i,j,v)\Q) = \<i>x,i,j,v)\y,{H idfV ) Xt y). 

We may obtain the unitary U g in the usual way by applying U g , copying the output, then applying C/J poj . 

Using the method of Ref. 0, the Hamiltonian H^ j ^ may be simulated using a call to U g and a call to U^. As 
z n is of order log* n, the number of black-box calls to / for the simulation of each Huj v -\ is 0(log* n). Using these 
values, along with Eq. , we obtain the number of black-box queries as in Eq. © • □ 

Another issue is the number of auxiliary operations, which is the number of operations that are required due to the 
overhead in calculating T(x,i, j). It is necessary to perform bit comparisons between a maximum of z n + 2 numbers 
in the first step, and each has n bits. This requires 0(n log* n) operations. In the next step the number of bits is 
0(log 2 n) bits, which does not change the scaling. Hence the number of auxiliary operations is 



O (n(log* n) 2 d 2 5 2fe (dV) 1+1/2fc A 1/2fc 
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This scaling is superior to the scaling n 10 in Ref. 

Next we consider the error introduced by calculating the matrix elements to finite precision. Given that the matrix 
elements are represented by 2n' bit integers, the error cannot exceed ||/2™ . The error in calculating exp(— iHnj jV -\t) 
will not exceed r/2™ so the error in the integrator due to the finite precision does not exceed 4m5 fe r/2 n . This 
error can then be kept below e/2 by choosing 

n' > 5 + log 2 (rd 2 5 fc /e). 

The total error may be kept below e by choosing the integrator such that the integration error does not exceed e/2. 



VI. CONCLUSIONS 

We have presented a scheme for simulating sparse Hamiltonians that improves upon earlier methods in two main 
ways. First, we have examined the use of higher order integrators to reduce the scaling to be close to linear in ||.ff||f. 
Second, we have significantly improved the algorithm for the decomposition of the Hamiltonian, so the scaling of the 
number of black-box calls is close to log* n, rather than polynomial in n. In addition we have shown that the scaling 
cannot be sublinear in \\H\\t (for reasonable values of n). 
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